function [P, F, dfR, dfF, x1, yp1, x2, yp2, xf, yf] = bears_f_test(x1r, y1r, x2r, y2r)
% Use F test to distinguish 2 models, either
% the data x1r, y1r and x2r, y2r is fit best by a
% common line, or by 2 separate lines
% Return the values P, F and df.
% after Bear et al., 1992.
%
% clean up the input data sets first. Polyfit doesn't like NaN's in the inputs
x1=x1r(find(~isnan(x1r) & ~isnan(y1r)));
x2=x2r(find(~isnan(x2r) & ~isnan(y2r)));
y1=y1r(find(~isnan(y1r) & ~isnan(x1r)));
y2=y2r(find(~isnan(y2r) & ~isnan(x2r)));

% now fit line to first set

[p1, s1] = polyfit(x1, y1, 1);
yp1 = polyval(p1, x1);
SSE_1 = sum((y1-yp1).^2);
%
% fit line to second set
[p2, s2] = polyfit(x2, y2, 1);
yp2 = polyval(p2, x2);
SSE_2 = sum((y2-yp2).^2);
%
SSE_F = SSE_1+SSE_2; % full model
%
% fit line to both (reduced model)
[pf, sf] = polyfit([x1 x2], [y1 y2], 1);
yf = polyval(pf, [x1 x2]);
SSE_R = sum(([y1 y2] - yf).^2);
%
% now compute the degrees of freedom for the 2 models
%
dfF = length([x1 x2])-4;
dfR = length(x2)-2;
F = abs(((SSE_R - SSE_F)/(dfR-dfF))/(SSE_F/dfF)); % returns F value for models.
P = fdist(F,2,dfF); % get the P value associated with the F value

% for convience, we sort the big result here -= can use for plotting later.
[xf, ix] = sort([x1 x2]);
yf = yf(ix);

return;
